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Abstract 

A two-dimensional Navier-Stokes solver developed for 
detailed study of wave rotor flow dynamics is described. 
The CFD model is helping characterize important loss 
mechanisms within the wave rotor. The wave rotor 
stationary ports and the moving rotor passages are 
resolved on multiple computational grid blocks. The 
finite-volume form of the thin-layer Navier-Stokes 
equations with laminar viscosity are integrated in time 
using a four-stage Runge-Kutta scheme. Roe’s 
approximate Riemann solution scheme or the 
computationally less expensive Advection Upstream 
Splitting Method (AUSM) flux-splitting scheme is used to 
effect upwind-differencing of the inviscid flux terms, 
using cell interface primitive variables set by MUSCL- 
type interpolation. The diffusion terms are central- 
differenced. The solver is validated using a steady 
shock/laminar boundary layer interaction problem and an 
unsteady, inviscid wave rotor passage gradual opening 
problem. A model inlet port/passage charging problem 
is simulated and key features of the unsteady wave rotor 
flow field are identified. Lastly, the medium pressure 
inlet port and high pressure outlet port portion of the 
NASA Lewis Research Center experimental divider cycle 
is simulated and computed results are compared with 
experimental measurements. The model accurately 
predicts the wave timing within the rotor passage and the 
distribution of flow variables in the stationary inlet port 
region. 

Introduction 

The wave rotor is an internal flow device used to 
transfer energy between gas streams of differing energy 
density. The energy exchange is accomplished through a 
timed sequence of unsteady compression and expansion 
wave processes that take place within the wave rotor 
passages which are arranged circumferentially at the 
periphery of the cylindrical rotor as shown in Figure 1. 
The waves are initiated as the aids of the wave rotor 
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passages sequentially open (or close) to stationary inlet 
and outlet ports. At an inlet port the low pressure gas in 
a rotor passage is exposed to the high pressure "driver* 
gas and a compression wave moves into the passage, 
compressing the "driven" gas on-board the rotor passage. 
The contact discontinuity between the driver gas and the 
driven gas follows the compression wave as the driver 
gas enters the passage and an expansion wave propagates 
into the inlet port region. At an outlet port, the typically 
higher pressure passage gas is exposed to the low 
pressure port flow and an expansion wave moves into the 
passage, reducing the on-board pressure, and exhausting 
the on-board gas to the outlet port. The charging and 
exhausting processes are used in the wave rotor to effect 
different gas cycles, including three-port divider cycles 1 , 
like the NASA LeRC wave rotor experiment described 
below, in which a medium stagnation pressure stream is 
split into a high stagnation pressure stream and a low 
stagnation pressure stream, and four-port Brayton 
topping-cycles, with application in propulsion systems. 2 

The wave rotor presents a challenging application for 
high resolution Navier-Stokes schemes. The numerics 
must resolve and accurately track the shock waves and 
contact discontinuities and model strong vortical flows 
resulting from unsteady shock/boundary layer interaction 
and shock/contact interaction. A two-dimensional (0,z) 
thin-layer Navier-Stokes solver has been developed for 
detailed study of the wave rotor flow dynamics. The 
multi-block solver uses high resolution upwind- 
differencing of the inviscid flux terms using either the 
Roe approximate Riemann solution scheme 3 or the 
AUSM flux splitting scheme 4,5 and central-differencing of 
the diffusion terms. 

The present modeling effort supports fundamental 
experimental 6 and computational 7 * 8 wave rotor research 
which is currently underway at the NASA Lewis 
Research Center aimed at characterizing the loss 
mechanisms due to friction, leakage, and gradual passage 
opening-time effects. A detailed study of the effects of 
geometry (e.g., inlet port lip rounding and passage 
profiling) on the shape and evolution of the contact 
discontinuity, on vortical flow structures generated within 
the rotating rotor passages that lead to mixing losses, and 
on the unsteadiness and nonuniformity of the flow in the 
stationary ports was conducted using the solver described 
here. 9 
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A detailed overview of the CFD modeling of wave 
rotor flow dynamics in the (0,z)-dimensions is presented 
in this paper. A model wave rotor problem, specifically 
developed to study the inlet port/passage charging 
process is first described. ITiis description identifies 
some of the important features of the wave rotor flow 
field. The numerical modeling is then provided, 
followed by examples of CFD and experimental model 
validation. Computed results for the model inlet port 
problem are presented. The results have helped identify 
the key features of the wave rotor flow field. Finally, 
computed results of a simulation of the medium and high 
pressure port portion of the NASA LeRC divider cycle 
are compared with experimental measurements. 

Model Inlet Port Problem Description 

Figure 2 provides a schematic diagram showing a 
group of rotor passages moving past an inlet port. In 
this view (looking radially-inward), the wave rotor 
passages are unwrapped from the periphery of the rotor. 
The right boundary of the rotor passages is a wall. 

(Note that the portion of the NASA LeRC divider cycle 
considered later in this paper is similar to the model 
problem shown here; however, in the divider cycle, a 
portion of the right-hand wall is replaced by the outlet 
port.) 

During the charging process shown, a lower pressure 
"virgin" passage is exposed to a higher total pressure 
inlet port flow. As the virgin passage gradually opens to 
the higher total pressure inlet port flow, a compression 
wave, eventually steepening into a shock, forms at the 
leading end of the passage and moves into the rotor 
passage. This wave is followed by a contact 
discontinuity between the driver fluid and the driven fluid 
as shown in Figure 2. During the gradual opening, an 
expansion wave is generated and moves into the inlet 
port region, reducing the total pressure of the driver fluid 
as it enters the passage. The rightward running shock 
wave compresses the on-board (driven) gas and 
eventually reflects off the right-hand wall. The leftward 
running reflected shock recompresses the driven gas, 
interacts with the contact discontinuity, and compresses 
the driver gas within the passage. The inlet port is 
designed so that the rotating passages close just as the 
reflected shock reaches the left end of the rotor passage. 

Computational Model 

The unsteady quasi-3-D rotor viscous code (RVCQ3D) 
of Chima 10 provided a skeletal structure for the finite- 
volume code (WROTOR2) developed in this work. Both 
single-passage and multi-passage/inlet port computations 
have been used to model the wave rotor flow field. In a 
single-passage calculation, the flow in a single rotor 
passage is modeled as that passage moves past either 
walls or inlet and outlet ports. 


An example multi-block computational domain is 
shown in Figure 2 and includes the inlet port region and 
the multiple rotor passages that can influence flow in the 
inlet port. The inlet port is modeled with a stationary 
computational block while the moving rotor passages are 
each resolved on separate, though identical, 
computational blocks. The walls separating the 
individual passages are infinitely dun in the current work 
except at the leading edge where profiling the passage 
wall can be included. To date, any outlet ports (e.g., 
like those of the NASA LeRC divider cycle discussed 
below) influencing a simulation are simply represented by 
uniform static pressure surfaces rather than by separate 
computational blocks. Modifying the multi-block solver 
to include outlet port blocks is a straight-forward 
extension of the work presented here. 

Model Restrictions 

In general, centrifugal forces acting on the gas within 
the rotor passages work to establish nonuniform, radially 
stratified flow. The centrifugal effects provide a 
mechanism for distortion or instability (and mixing) at 
the contact discontinuity for driver and driven gases of 
differing density 11 ; however, the degree to which this 
mechanism is important might be expected to scale with 
the square of rotor speed (or wheel Mach number) and 
the relative height of the rotor passage to the radius of 
the rotor. Influenced by the relatively low rotor speeds 
(4000 to 8000 rpm) and the dimensions (Ar/R — 0.067) 
of the NASA LeRC experiment, in this study both 
centrifugal and Coriolis forces are neglected by choosing 
to model the rotor passages and inlet ports in the 0-z (or 
y-x) coordinates shown in Figure 2. Leakage flows, both 
from the rotor passages to the space between the rotor 
and the outer casing, and from passage to passage, are 
not yet included in this 2-D model. Fluid-to-wall heat 
transfer is another potential loss mechanism neglected in 
the calculations presented here as the blades separating 
the passages are treated as adiabatic surfaces. 

Governing Equations 

The unsteady, compressible 2-D thin-layer Navier- 
Stokes equations with laminar viscosity are used to model 
the wave rotor working fluid, which is treated as a 
calorically and thermally perfect gas. The governing 
equations in boundary-fitted coordinates (£, 77 ) under the 
thin-layer approximation 12 are 

dQ Be da 

dt dt, d 11 di\ 

where Q = (p, pu, pv, pe) T /J is the vector of conserved 
variables, F and G are the inviscid flux vectors given by 
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U = V{* (u-i) and V - Vfj* (ji-i) are the £ and 7 ? 
normalized contravariant velocity components relative to 
the grid motion, u and v are the axial and tangential (x 
and y) components of velocity, z is the grid velocity, p 
is the fluid mass density, p is the static pressure, e is the 
total specific energy, h is the total specific enthalpy, 

(i'JU ) x + t^v , and = x^u + x^v. The 

shear stress (e.g., r^) and heat flux components (e.g., 
q") are also expanded in the boundary-fitted coordinates 
and, under the thin-layer approximation, all derivatives 
with respect to the streamwise coordinate ({) are 
neglected. The normalized metric terms (e.g., 

% x * £ x / |V$ |) and the surface areas of the finite- 
volumes (e.g., |V{ | / J where J is the Jacobian of the 
(x,y) ** ( $ / 11 ) transformation) are calculated using 
the specified grid. In this work, the grid speed metric { t 
is zero while t| c = -r<o where r is the radius at which 
the calculations take place (e.g., the passage average 
radius) and 0 ) is the angular velocity of the rotor. 

Time Integration 

The governing equations in finite-volume form are 
integrated in time using a second-order accurate explicit 
four-stage Runge-Kutta scheme (see Jameson et al. 13 ). 
The time step is dictated in part by the CFL number, a, 
which is nominally set to 1.2. For solutions of steady 
flows or for unsteady inviscid flow fields resolved on 
nearly uniform grids (on which the CFL number is 
nearly uniform in space), the implicit residual 
smoothing 14 

( 1 (1 -€ n 5 nn ) R k =R k 

where 5^ and 5^ are central-difference second 
derivative operators, can be applied after each k* Runge- 
Kutta step to smooth the residuals R k ~ (dQ/dt) k so 
that the local CFL number can be increased. In this 
work the smoothing coefficient suggested by Turkel and 
van Leer 15 , e 5 = 0,9o/8 is used for the {-direction 
smoothing and * 0 . 5 is used for the 77 -direction 
smoothing. Typically, o = 3 is used when the residuals 
are smoothed. The residual smoothing adds 
approximately a 3 to 5% overhead to the overall cpu- 
time. 


Spatial Discretization 

In the cell-centered finite-volume formulation, the 
inviscid and viscous flux vectors are calculated at the cell 
interfaces which are constructed using lines connecting 
the grid points. The shear stress and heat flux 
components of the viscous flux vector are central- 
differenced at the cell interfaces using the primitive 
variables of the adjacent cells and transport properties 
based on the arithmetic mean of the static temperatures in 
the adjacent cells. The inviscid flux vectors are 
computed at the cell interfaces so as to effect a second- 
order accurate upwind-difference using the MUSCL-type 
interpolation and difference schemes described below. 

Interpolation. Spekreijse’s interpolation 16 of the 
primitive variables (p,u, p) with Van Albada limiting 17 is 
used to set the "left" and "right" states of the cell 
interfaces. Our choice of MUSCL-type 18 interpolation to 
effect a second-order accurate, monotonicity preserving 
scheme, was influenced by the work of Anderson et al. 19 
As implemented by Orkwis and McRae 20 , the primitive 
variables are set (e.g., at a constant -77 interface) by 

1 T 1 (2db 2 + 8 2 ) da+ ( da 2 + 2 6 2 ) db 

i+ *l i 2 2db 2 - dadb + 2 da 2 +3 5 2 

R _ 1 (2db 2 + 6 2 ) dc+ (dc 2 + 2 8 2 ) db 

i4 i ~ i+1 2 2dfo 2 - dbdc + 2 dc 2 + 3 6 2 


where da = , db = u i+1 - u d , and 

dc = u i+2 -u i+1 , 5 = 1 x 10 -5 , and where "u" here 
represents the primitive variables. Having set the left 
and right states, the inviscid flux at a cell face is 
calculated by using the approximate Riemann solution of 
Roe 3 or the flux-splitting scheme of Liou and Steffen 4 . 

Roe flux difference splitting. The approximate 
Riemann solution of Roe 3 provides an estimate of the 
inviscid flux at a cell face of constant -77 (for example) 
through 

= -|<Z(.U £ ) +E(u r ))- I |||8£ 

= ■§<**+*; - it 1 

^ is 1 


where A denotes Roe-averaged variables 3 , "L" and "R" 
denote the left and right states of the Riemann problem 
set by the interpolation, and E z and E R are the flux 
vectors computed using the primitive variables (jJl and 
Ur) of the left and right states. is the 4 modified 
eigenvalue and t 1 is the P right eigenvector of the 
Jacobian matrix A = dE/dQ (see, for example, Hirsch et 
al. 21 ), constructed using Roe-averaged variables, and &w 1 
is the jump (left to right state) in the l 111 characteristic 
variable. The modified eigenvalues are the eigenvalues 
of the Jacobian £-i.e., ^ = X 2 = &|V$ | /J and 



vt 3 , 4 - ( U± £) |V£ \/J (see Steger 22 and his reference to 
Warming et al. 23 for the eigenvalues in moving 
generalized coordinates)— corrected using the entropy 
condition fix of Harten and Hyman 24 to assure that 
expansion shocks are not admitted by the solution. The 
generalized coordinate formulation of Roe’s flux 
difference splitting used in this work is that presented by 
Grossman and Walters 25 and Walters and Thomas. 26 The 
flux vector in the | -direction (at a cell face of constant-^) 
can be written as 

where 


AUSM flux vector splitting. As an alternative to the 
computationally expensive Roe scheme, the Advection 
Upstream Splitting Method (AUSM) of Liou and Steffen 4 
(see also Liou 5 ) has also been investigated to set the 
inviscid flux vector. Accommodating the grid motion, 
and in generalized coordinates, the AUSM splitting is 
used in this work by calculating the inviscid flux (e.g., 
here at a cell face composed of a constant-77 line and 
using the notation of Liou 5 ) as 

-^{i - 1 *-*4} 

where "L" and "R" again denote the left and right states 
respectively and where 
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where u 5 = U+V%-x = U-% t , the jumps are from the left 
to right states (e.g., 6 u = u R -u L ), and where the 
normalized metrics are those of the cell surface, 
described earlier. Note that 8u* = bU for the case 
where the grid velocity x is uniform in space, as is the 
case in the wave rotor. The G-vector (e.g., is 

computed at the cell faces composed of constant-! lines 
in an analogous fashion. 

We note that the grid motion influences the above 
formulation in that the contravariant velocity components 
(U and V) are those relative to the moving grid and 
therefore the eigenvalues are affected by the grid motion. 
Because the Jacobian A is given by A = € t Z * AL where 
AL is the Jacobian dE/dQ for a stationary grid system, 
the eigenvectors of the Jacobian AL of the stationary 
system are also the eigenvectors of the Jacobian A of the 
moving grid system; i.e., the Roe dissipation term 
(.1 |^|5j2) on the moving grid system can be decomposed 
into the standard (fixed grid system) right eigenvectors 
and jumps in characteristic variables, as done above, and 
the eigenvalues are shifted by the grid speed. 


given that Ui=Ul + U R and Pi=Pl + Pr. Defining 

the local contravariant Mach number of a given state as 

8- &/a, where a is the local speed of sound, the Ui 

and pi are given by 
2 
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The AUSM scheme effectively sets the convective part of 

the inviscid flux at the cell surface to be U± $£ if 

2 

Ux ^ 0 or Ui otherwise. Liou and Steffen 4 note that 
2 "2 

the third-order pressure splitting used in the subsonic 

region can cause small oscillations immediately (one cell) 

downstream of a shock. A simpler first order expansion 

can be used in the subsonic region 4 : p ± = p ( 1 ± Af) /2 . 

The third-order pressure splitting is used in the present 

work. Liou and Steffen also point out that the AUSM 

scheme, like van Leer’s flux splitting scheme, uses a 

scalar coefficient (Ux) and therefore requires only 0(n) 

2 

operations for an n-dimensional flux vector as compared 
to 0(n 2 ) operations in the Roe scheme because of the 
matrix coefficient. For example in a computation in this 
present work of Mach 2 steady inviscid flow over a 10% 
bump, 1000 iterations with the AUSM scheme required 
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0.70 time units as compared to the 1.0 time units 
required by the Roe scheme with the entropy fix. The 
AUSM scheme shows promise of computational 
efficiency while maintaining accuracy. 

Hie overall thin-layer Navier-Stokes solver is spatially 
and temporally second-order accurate. It has been 
demonstrated for both the Roe flux difference splitting 
(see, for example, van Leer et al. 27 and von Lavante 28 ) 
and the AUSM splitting scheme (see Liou 5 )that the 
combination of the upwind treatment of the inviscid flux 
terms and the central-difference of the diffusion terms 
accurately models steady viscous flow fields. Here we 
apply these methods to the unsteady wave rotor viscous 
flows. 

Boundary and Interface Conditions 

The unsteady wave rotor flow field has been modeled 
using both a single-passage solver and a multi-block 
solver. 

Single-passage. In the single-passage computations, 
the flow field is computed in the frame of reference of 
the moving passage and moving boundary conditions are 
applied at the ends of the single computational block as 
shown schematically in Figure 3. In the single-passage 
computations reported here, the inflow boundary is set by 
completely specifying the primitive variables at the 
inflow plane (State a) and by requiring that the expansion 
fan emitted upstream, from the passage, be infinitely 
weak (cf. Eidelman 29 ); i.e., the primitive variables at the 
phantom point (i = 0) are those of State a. The v 
component of velocity is set to the rotor speed. At an 
outflow boundary (at the right end of the passage shown 
in Figure 3), the static pressure is specified at the 
outflow plane (state b) and the invariance of 
R+ = 2a/ (y-1) + u and p/p y across the fan is used to 
set the phantom point (implying that the emitted contact 
and shock are infinitely weak), and v is linearly 
extrapolated from the interior. The effects of the gradual 
opening of the rotor passages to the inflow or outflow 
ports is simulated, after Eidelman 29 , by gradually 
replacing the wall boundary conditions described next by 
the appropriate inflow or outflow conditions, based on 
the rotor linear speed. 

At constant-^ walls, dp/dn = 0 , where n is the local 
normal to the wall, and V = 0 is enforced within the 
inviscid flux vectors. At constant-^ walls, dp/ dn = 0 is 
enforced by using phantom points in the wall; e.g., 

U 0 = -U 1 and p 0 = p 1 are required and additionally, 

V Q = V ± is required for an inviscid flow. For viscous 
calculations, dT/dn = 0 and U=V= 0 are also required 
at all the walls and these conditions are enforced by 
using phantom points. 

Multi-block. In the multi-block calculations, the 


phantom points to the left of the inlet port inflow plane 
are set by specifying the total pressure and total 
temperature and by extrapolating the left-running 
Riemann invariant R~ * 2a/ (y-1) ~ ^ (where is 
the covariant ^-velocity component) from the first 
interior point. The downstream boundary of the inlet 
port and the upstream boundary of the rotor passages 
currently within the span of the inlet port coincide as 
shown in Figure 2. This surface is treated as an 
interface and is updated at the end of a complete time 
step, after incrementally advancing the grid as 
appropriate for the current time step and the rotor linear 
speed. The phantom point at the downstream boundary 
of the inlet port is set using linear interpolation of the 
conserved variables on the current rotor solution while 
the phantom point at the inlet of the rotor passages is set 
using linear interpolation of the conserved variables on 
the current inlet port solution. Outside of the inlet port 
span, the phantom points at the left end of the rotor 
passages are set to effect the appropriate wall boundary 
conditions as described above for the single-passage 
calculations. 

At the right-hand end of a rotor passage, the phantom 
points are set to effect the constant-1 wall boundary 
conditions for those points outside of an outlet port span. 
For simulations in which outlet ports are included, the 
phantom points exposed to the outlet port flow are set as 
in the single passage calculations; i.e., a uniform static 
pressure is specified at the outlet boundary and the 
invariance of the R+ and p/p 1 across the fan that is 
assumed to enter the passage is used to set the phantom 
point, and the v velocity component is linearly 
extrapolated from the interior; however, if the resulting u 
velocity component at the phantom point is negative, 
indicating local inflow rather than outflow, the phantom 
point stagnation temperature is then specified, the u 
velocity component is set to zero, and the mass density is 
calculated based on the local temperature and the same 
specified back pressure. 

Model Validation 

Steady Flow 

After MacCormack and Lomax 30 , the ability of the 
solver to model accurately a stationary oblique 
shock/laminar boundary layer interaction is evaluated by 
comparing computed results of the present work with the 
experimental data of Hakkinen et al. 31 and calculated 
results of Liou and Steffen 4 using their AUSM scheme. 
Figures 4a and 4b show comparison of the normalized 
static pressure (p val2 / P~, in) and the skin friction 
coefficient ( C t s 2x v / (pu 2 )„ fin ) of a Mach 2 flow 
along a flat plate (with a leading edge at x = 1). Both 
the AUSM flux-splitting and Roe flux-differencing 
schemes as implemented in the present work are 
compared. The present computations were performed on 
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a 121 x 101 stretched grid. Liou and Steffen resolved 
the flow field with a 75 x 65 stretched grid. A 
32-degree oblique shock interacts with the laminar 
boundary layer at x=2 where the Reynolds number is 
2.96 x 10 s , and causes separation. The computed results 
using the Roe and AUSM schemes as implemented in the 
present study are virtually identical. The static pressure 
distributions predicted by both the Roe and AUSM 
schemes in this present study agree very well with the 
experimental data of Hakkinen et al. 31 The static 
pressure predicted by Liou and Steffen using their 
AUSM scheme is low in the recirculation region (as 
discussed by Liou and Steffen 4 ). The skin friction 
coefficient predicted by the present method is lower than 
the experimental data over much of the plate. The 
results of Liou and Steffen’s calculation better predict the 
extent of the recirculation zone indicated by the skin 
friction coefficient. Qualitatively, the. CFD predictions of 
the present study and experimental data of Hakkinen et 
al. match well. 

Unsteady Flow 

The wave rotor passage gradual opening problems 
considered by Eidelman 29 have been repeated as a form 
of CFD validation. Eidelman used a second-order 
extension of the Godunov method (see Eidelman et al. 32 
for details of their model). Figure 5 provides a 
comparison of snapshots in time of the static pressure 
contours reported by Eidelman 29 and computed here as a 
single passage calculation. The problem models a single 
wave rotor passage moving downward in Figure 5 at 200 
m/s. The present calculations use the AUSM flux- 
splitting scheme and are performed on a uniform grid of 
121 x 25 points which provided the same resolution as 
that used by Eidelman. The conditions in the inlet 
port/passage interface are p b = 0.182 MPa (1.80 atm), 

11 * = 150 m/s, v* = 0, and p m = 1.81 kg/m 3 , and the 
conditions on board the passage prior to opening are p> 

= 0.103 MPa (1.0 atm), p© = 1.2 kg/m 3 , and u© = 0. 
These conditions assure that there is an insignificant 
density change across the driver gas/driven gas contact. 
The passage itself does not move in Eidelman’ s problems 
but rather the gradual opening of the passage to the inlet 
port flow is simulated by gradually replacing the left wall 
boundary condition with the inlet port flow conditions in 
accordance with the rotor speed of 200 m/s such that the 
passage is fully open at 0.1 ms. The time sequence 
compared in Figure 5 indicates that the two models agree 
well. The flow field appears to better resolved in the 
current study. Both models clearly show the 
development of the compression wave from an initial 
"spherical" wave, through reflections of the pressure 
wave first off the trailing (upper) passage wall (evident at 
t=0.086 ms), then the off the leading wall (evident at 
t=0.127 ms and t= 0.168 ms), and finally off the trailing 
wall again (evident at t=0.252 ms). A discussion of the 
reflection of the compression wave and the highly 
rotational flow that results from the gradual, rather than 


instantaneous passage opening can be found in 
Eidelman’s paper. 29 

The single-passage problem was extended to model the 
shock/contact interaction and shock/boundary layer 
interactions by including viscous effects and by allowing 
the initially rightward running shock wave to reflect off 
the right-hand wall. The inlet port/passage interface 
conditions in this computation are p in = 0.289 MPa (2.85 
atm), u* = 283 m/s, and p m = 4 kg/m 3 , and the rotor 
speed in this calculation is 250 m/s. The on-board 
passage conditions are the same as in the inviscid case. 
These conditions are those of a second case considered 
by Eidelman 29 and assure that the driver gas/driven gas 
contact discontinuity is significant so that the 
shock/contact interaction can be studied. The 
calculations were performed on a 121 x 61 stretched grid 
using the AUSM flux-splitting. Figure 6 shows velocity 
vectors and density contours at a time when the reflected 
shock, apparent at the left side of the picture, has moved 
through the contact discontinuity. (Note that in this 
picture the left wall was gradually withdrawn in the 
direction opposite of that in the inviscid problem above.) 
The strong vortical structure following the shock/contact 
interaction is clearly visible. Note that the highly 
rotational flow leads to mixing between the driver and 
driven gases and is a loss mechanism in the wave rotor. 
At the position of the reflected shock, a strong unsteady 
shock/boundary layer interaction is evident. The 
boundary layer fluid into which the shock moves lacks 
the momentum to negotiate the shock pressure rise so 
that separation occurs, a recirculation region is formed, 
and the upstream boundary layer is lifted, resulting in 
turning of the (relative) supersonic free stream through 
an oblique shock stem. The flow is then turned again to 
the axial direction via another oblique shock stem. The 
reflected shock is therefore a running lambda-shock. 

From Figure 6 it is apparent that the effective distance in 
which the upstream fluid is brought to rest is the length 
of the recirculation region. 

Model Inlet Port Simulation 

Multi-block calculations have been used to study the 
influence of inlet port comer rounding and of profiling 
the leading passage wall on the shape and behavior of the 
contact discontinuity established by the charging process 
at an inlet port. More detailed results of the study have 
been reported in Ref. 9. A brief description of the 
model inlet port simulation is presented below along with 
example results showing the evolution of the contact 
discontinuity during the rotor charging process. 

The inlet port opening of this simulation spans five 
rotor passages, each of which has a blade-to-blade width 
of 0. 1 units and a length of 1.0 units. The inlet port 
walls are each slanted at 30-degrees (off the horizontal) 
so that the mean y-component of velocity nearly matches 
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the rotor linear speed. Seven rotor passages are included 
in a calculation at any given time. As a virgin passage 
enters the computational domain, the most advanced 
passage is discarded. At the inlet port inflow plane, the 
stagnation pressure and temperature are specified as 3.0 
and 1.0, respectively, in nondimensional units. These 
conditions are kept constant in space and time. The on- 
board virgin passage stagnation pressure is 1.0, the on- 
board stagnation temperature is 1.0, and the relative 
velocity is zero. The wheel speed is selected so that the 
initially rightward running shock wave in a given passage 
reflects off the right-hand wall and arrives back at the 
inlet port/passage interface plane as that passage is 
closing. The right-hand side of the rotor passages move 
past a fixed wall everywhere in the example simulation. 
The inlet port grid is 56 x 101 points and each of the 
seven passage grids has 106 x 31 points, for the inviscid 
flow shown here. The AUSM flux-splitting scheme was 
used in these computations. 

The initial conditions were uniform Mach number (0.5) 
in the inlet port and virgin passage conditions in each 
passage. The inlet stagnation conditions and virgin 
passage conditions are constant and therefore the 
unsteady flow field reaches temporal periodicity after a 
large number of rotor passages have traversed the inlet 
port. The unsteady flow field typically begins to assume 
periodic behavior in about three times the time it takes a 
given passage to traverse the inlet port span (e.g., three 
thousand iterations). True periodicity however requires 
longer r unnin g time and is dependent on the grid 
resolution. Ten to fifteen thousand iterations were 
required to reach periodicity in these inviscid flow 
simulations. 

Sharp Inlet Port Lip 

Figure 7 shows density contours for an inlet port with 
a sharp cornered lip (no rounding). The passage is 66 
percent open at the time shown. A rightward running 
compression wave in the first (bottom) passage is seen to 
steepen quickly into a rightward running shock wave (as 
seen in the second passage). Simultaneously an 
expansion wave enters the inlet port region. The 
rightward running shock wave compresses the on-board 
driven gas, setting it into motion and allowing the driver 
gas to enter the passage. The contact discontinuity 
moves with the local fluid velocity. Due to the sharp 
comer inlet lip and the asymmetry of the gradual opening 
process, the driver fluid is initially (see the first passage 
at the bottom of Figure 7) concentrated along the leading 
wall of the passage. Mach contours (not shown here) 
show that the driver fluid behind the most forward part 
of the interface is supersonic and is running faster than 
the subsonic driven gas by which it travels. There is 
large vorticity at the contact, apparent in the subsequent 
tumbling of die interface as seen in the first, second, and 
third passages. Note that a part of the driven gas is 
"captured" by the tumbling interface. 


The initially rightward running shock reflects off the 
right wall and recompresses the driven gas, bringing the 
driven gas to rest (see the third passage). The leftward 
running reflected shock then interacts with the contact 
discontinuity as seen in the fourth passage. Note that any 
protuberance in the contact is exaggerated by the 
shock/contact interaction. The reflected shock further 
compresses the driver fluid (see the fifth passage). The 
reflected shock returns to the inlet plane and is emitted 
from the passage into the inlet port as the passage closes. 
The emitted shock wave is visible in the inlet port region 
shown in Figure 7. This fluid leaving the sixth passage 
is accelerated around the passage wall and enters the fifth 
passage. The flow in the inlet port is unsteady due to the 
emitted shock wave and the leftward moving expansion 
wave created as new passages open to the inlet port. 

Inlet Lip Rounding 

Figure 8 shows the density contour plots for a wave 
rotor inlet port with lip rounding at the same time as the 
previous example. Comparing Figures 7 and 8 shows 
that the comer rounding has a large effect on the 
behavior of the contact discontinuity. In this case the 
leading lip of the contact discontinuity is located at the 
trailing wall of the passage (see the first, second, and 
third passages). The contact discontinuity curls in the 
direction opposite that of the sharp lip case after the 
reflected shock/contact interaction. Due to the 
asymmetry of the gradual passage opening process, a 
portion of the driven fluid is again captured by the driver 
fluid; however, the comer rounding appears to have 
significantly reduced the amount of driven fluid captured, 
suggesting that comer rounding can help control 
driver/driven gas mixing. 

NASA LeRC Divider Cycle Simulation 

The medium pressure inlet port and the high pressure 
outlet port region of the NASA LeRC three-port divider 
cycle experiment (see Wilson and Fronek 6 ) has been 
modeled. In a three-port divider cycle, typically a small 
fraction (e.g., 0.3) of the medium pressure stream that 
enters through the inlet port is compressed to a higher 
stagnation pressure (e.g, P high / Pmeditm ~ L 3) and is 
then exhausted into the high pressure outlet port. The 
remainder of the stream is expanded later in the cycle to 
a lower stagnation pressure (e.g., p low / P medium ~ 0.1%) 
and is exhausted to the low pressure outlet port (see, for 
example, Kentfield 1 ). In this simulation, the medium 
pressure port and 39 rotor passages that are near the 
medium pressure inlet port or the high pressure outlet 
port are computed with individual blocks. This portion 
of the NASA LeRC divider cycle is similar to the model 
inlet port problem described earlier (and shown 
schematically in Figure 2); however, the high pressure 
outlet port on the right end of the rotor allows the twice- 
compressed driven flow to leave the rotor passage. The 
outlet port is represented here by a uniform static 
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pressure surface. A more detailed discussion of this 
simulation has been reported in Ref. 9. 

Two viscous calculations for the divider cycle were 
performed using Roe’s scheme. In the first, the virgin 
passage flow is set to zero relative velocity and the on- 
board pressure and temperature are set to best estimate 
conditions based on experimental data. In the second 
calculation, the virgin passage conditions are a 2-D 
construction of the passage conditions just prior to 
opening to the inlet port based on 1-D computations of 
the full divider cycle. 8 

Experimentally measured quantities (private 
communication with J. Wilson at NASA LeRC, 1993) 
are specified at the inlet and exit ports. The inlet port lip 
of the NASA experiment is rounded in order to prevent 
flow separation as the fluid enters the moving virgin 
passage. The inlet port spans 27 passages at a given 
time. The walls of the inlet port are parallel and are 
oriented at 22-degrees off the horizontal. A 20.3 cm 
(8.0 in.) long by 17.1 cm (6.75 in.) wide (span) portion 
of the inlet port is resolved by a 81 x 101 stretched grid. 
The rotor passages are 0.64 cm (0.25 in.) wide by 45.7 
cm (18.0 in.) long and each of the 39 rotor passages is 
resolved by a 151 x 15 stretched grid. The average 
passage radius is 15.2 cm (6.0 in.) and the rotor spins at 
3924 rpm. 

Figure 9 shows static density contours at a point when 
the first passage is 66 percent open to the inlet port. The 
calculation is performed using the 1-D CFD-based virgin 
passage initial conditions mentioned above. The initial 
rightward running compression shock is mistimed and 
arrives late at the right-hand end of the passage; i.e., the 
passage has already opened to the outlet port flow, flow 
reversal occurs, and fluid enters the passage from the 
outlet port until the shock arrives. The shock then 
reflects not off the wall but off of the fluid entering from 
the outlet port. The late arrival of the shock is 
experimentally verified by a pressure transducer trace. 
Behind the rightward running shock is a rightward 
moving region of low density fluid that entered the 
passage earlier by leakage from the rotor casing. This 
low density fluid is evident in this 2-D calculation only 
because it is included in the virgin passage density 
distribution. Following this low density fluid is the 
contact discontinuity. The reflected shock passes through 
the contact, resulting in stretching of the interface, 
consistent with the inviscid and viscous results of the 
previous section. The reflected shock is mistimed and 
arrives at the inlet port/passage interface well before the 
port closes, consistent with experimental measurements. 

A rightward running expansion wave is generated as the 
passages exposed to the inlet port close. Likewise, a 
weak hammer shock is created at the right-hand end as 
passages emptying into the outlet port close. The effect 
of the comer rounding is evident in the low density 


region just above the rounding. 

Figure 10 compares the effectively time-averaged 
experimental static pressure distribution (private 
communication with J. Wilson at NASA LeRC, 1993) 
along the inlet port circumference (here normalized to 
one) with results from the two computations mentioned 
above. Qualitatively the CFD and experimental 
measurements agree fairly well; however, in both CFD 
calculations the static pressure predictions are low, 
indicating that the flow is moving relatively faster in the 
CFD solution than in the experiment. Indeed Figure 11 
shows that the CFD predictions of the total velocity are 
significantly higher than the experimentally measured 
values. It is noted that the effect of blockage due to the 
thickness of the passage walls, comprising nearly 16 % of 
the total port flow area, is not accounted for in the 
calculations for which infinitely thin walls have been 
assumed. Accounting for this blockage should improve 
agreement. Figure 12 compares the CFD predictions and 
the experimental measurements of the flow angle (tan 
[ (v/u)) circumferentially along the inlet port. The 
predicted flow angles are consistently higher than the 
measured flow angles; however, the qualitative behavior 
is effectively modeled. Finally, note that the results of 
the two CFD predictions qualitatively agree, suggesting 
that the discrepancy between the flow fields is due 
largely to differences in the level of the on-board 
pressure in the virgin passage rather' than due to the band 
of leakage fluid. This observation gives confidence in 
the ability to accurately simulate the qualitative behavior 
of the medium and high pressure port flow dynamics 
using simple, best estimate, quiescent flow conditions in 
the virgin passage. The agreement between the 
experimental data and the CFD prediction of the nearly 
steady port flow and between the CFD -predicted wave 
timing and that of the experiment, give confidence that 
the solver is doing a good job predicting the wave rotor 
flow dynamics in the vicinity of the medium and high 
pressure ports. 

Summary 

A 2-D, multi-block, thin-layer Navier-Stokes solver, 
developed to simulate wave rotor flow dynamics, has 
been described. The governing equations in finite- 
volume form are integrated in time using a four-stage 
Runge-Kutta scheme. MUSCL-type interpolation sets the 
"left" and "right" states used by either the Roe's flux 
difference splitting scheme or the AUSM flux-splitting 
scheme to compute second-order accurate, monotonicity- 
preserving, upwind-differences of the inviscid flux 
vectors at the finite-volume cell interfaces. The diffusion 
terms are central-differenced. The 2-D solver, 
WROTOR2D, has been validated using steady and 
unsteady flow problems. The AUSM scheme is found to 
require 0.7 times the cpu-time of the Roe scheme with 
entropy fix, as implemented in this work. The AUSM 
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scheme promises both computational efficiency and 
accuracy. 

A comparison of computed results, using both the Roe 
and AUSM schemes, with experimental data indicate that 
the solver accurately models stationary shock/laminar 
boundary layer interaction. A wave rotor passage 
gradual opening time problem was simulated using the 
AUSM scheme and computed results compared well with 
those of a separate study 29 . It is reasonable to expect that 
the success in predicting unsteady inviscid flow and the 
steady shock/boundary layer interaction problem carries 
over to the moving lambda-shock/boundary layer 
interaction problem. The unsteady shock/boundary layer 
interaction and the shock/contact interaction are common 
features of the wave rotor passage flow field, and were 
well resolved by the single-passage computation. 

A model inlet port/rotor passage charging problem was 
simulated using the AUSM scheme. Key features of the 
unsteady flow field, including the development of the 
main shock wave and the contact discontinuity, the 
reflected shock/contact interaction, and the emission of 
the reflected shock from the rotor passage into the inlet 
port region, were well resolved. An example of the 
influence of inlet port lip rounding on the shape and 
evolution of the contact discontinuity was provided. The 
sources of flow unsteadiness in the inlet port region were 
identified from the inlet port/passage charging problem 
results. 

Lastly, the medium pressure inlet port and high 
pressure outlet port portion of the NASA LeRC divider 
cycle experiment was simulated using the Roe scheme. 
The viscous flow field was computed on forty 
computational blocks, representing the stationary inlet 
port and thirty-nine moving rotor passages. Comparisons 
between computed results and experimental data were 
made. The computed wave timing is qualitatively 
verified by experimental pressure traces. The computed 
and measured static pressure, total velocity, and flow 
angle distributions along the inlet port qualitatively agree 
well. 
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Figure 1. Schematic diagram of a four-port wave 
rotor (courtesy of D. E. Paxson 7 ). 



Figure 2. Model problem schematic diagram showing 
multiple wave rotor passages moving past inlet port 
Computational domain denoted by dashed lines. 



INFLOW PLANE OUTFLOW PLANE 


Figure 3. Single-passage inflow and outflow boundary 
condition diagram. 
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a. Normalized static pressure b. Skin friction coefficient 

Figure 4. Comparison of oblique shock/laminar boundary layer interaction measured and computed 
distributions of a.) normalized static pressure, and b.) skin friction coefficient along flat plate. 
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a. Eidelman’s results. 29 b. Present study results. 

Figure 5. Comparison of computed static pressure contours by a.) Eidelman 29 , and b.) the present study 
for wave rotor gradual opening problem. 
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Figure 6. Density contours and velocity vectors of single passage viscous flow calculation showing the 
interaction of reflected shock and boundary layer and vortices generated by reflected shock and contact 
discontinuity interaction. 




Figure 7. Computed density contours of the model inlet port/passage charging problem 
with sharp corner inlet port lip. 
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Figure 8. Computed density contours of the model inlet port/passage charging problem 
with rounded inlet port lip. 
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Figure 9. Computed density contours of the medium pressure inlet port and high pressure outlet port 
portion of the NASA LeRC divider cycle experiment. - 
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Figure 10. Comparison of measured and computed static pressure circumferential distribution in 
the inlet port of the NASA LeRC divider cycle experiment. 



Figure 11. Comparison of measured and computed total velocity circumferential distribution in 
the inlet port of the NASA LeRC divider cycle experiment. 



Figure 12. Comparison of measured and computed flow angle circumferential distribution in 
the inlet port of the NASA LeRC divider cycle experiment. 
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